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ABSTRACT 

A  time-accurate  finite  element  model  for  simulating  the  fully- 
coupled  dynamic  response  of  flexible  multibody  systems  and 
liquid  sloshing  in  tanks  is  presented.  The  semi-discrete 
combined  solid  and  fluid  equations  of  motions  are  integrated 
using  a  time-accurate  parallel  explicit  solver.  The  FE  model 
consists  of:  hexahedral,  beam,  and  truss  solid  elements;  rigid 
bodies;  joints;  actuators;  hexahedral  incompressible  fluid 
elements;  and  quadrilateral  fluid-solid  interface  elements.  The 
fluid  mesh  is  modeled  using  a  very  light  and  compliant  solid 
mesh  which  allows  the  fluid  mesh  to  move/deform  along  with 
the  tank  using  the  Arbitrary  Lagrangian-Eulerian  formulation. 
The  fluid's  free-surface  is  modeled  using  an  acceptor-donor 
volume-of-fluid  based  algorithm.  The  motion  of  the  solid  and 
fluid  is  referred  to  a  global  inertial  Cartesian  reference  frame.  A 
total  Lagrangian  deformation  description  is  used  for  the  solid 
elements.  The  penalty  technique  is  used  to  model  the  joints. 
Numerical  simulations  are  presented  for  a  half-filled  tank 
supported  by  linear  springs  mounted  on  a  test  fixture. 

1.  INTRODUCTION 

Many  practical  applications  include  a  flexible  multibody 
system  carrying  liquid  filled  tanks.  The  multibody  system  can 
be  a  ground  vehicle  (truck,  train,  or  car),  ship,  airplane 
(commercial  or  military  jet,  helicopter,  etc.),  or  a  space 
structure  (space  station  or  satellite).  The  tank  can  be  a  payload 
tank  or  a  liquid  fuel  tank.  Optimum  system  design  in  those 
applications  requires  an  accurate  computational  model  of  the 
system  response,  which  in  turn  requires  accurate  modeling  of 
the  following: 

•  Incompressible  fluid  flow  in  a  moving/deforming  container 
including  accurate  modeling  of  the  free-surface,  turbulence, 
and  viscous  effects. 

•  Flexible  multibody  system,  including:  flexible  bodies,  rigid 
bodies,  joints,  frictional  contact. 

•  Coupling  between  the  solid  and  the  fluid  at  the  fluid- 
structure  interface. 

The  fluid  flow  is  governed  by  the  incompressible  Navier- 
Stokes  equations.  Finite  volume  [1,  2,  3],  finite  element  [4,  5, 
6]  or  particle  [7]  discretization  techniques  have  been  used  to 


model  the  fluid  flow.  In  addition,  many  techniques  for 

modeling  fluid  flow  with  a  free-surface  in  a  moving/deforming 

container  have  been  developed  in  the  literature,  these  include: 

(a)  Fixed  Cartesian  fluid  mesh  with  volume-of-fluid  (VOF) 
free-surface  model.  The  fluid  domain  is  a  Cartesian  mesh. 
The  container  moves  inside  this  mesh.  Therefore,  a  cut-cell 
technique  must  be  used  in  order  to  find  where  the  boundary 
of  the  container  intersects  with  the  fluid  cells.  The  cut-cell 
surfaces  are  then  used  to  impose  wall  impenetrability  and 
adhesion  boundary  conditions  [1,  2].  A  VOF  technique  [8, 
9]  is  used  to  model  the  liquid  free-surface  where  each 
element  has  a  VOF  value  between  0  (for  empty  elements) 
and  1  (for  elements  completely  filled  with  fluid).  The  free 
surface  is  reconstructed  for  each  element  using  piecewise- 
linear  planar  segments  using  the  VOF  value  of  the  element 
along  with  the  VOF  values  of  neighboring  elements  (which 
are  used  to  determine  the  normal  to  the  planar  surface). 
Note  that  when  VOF  algorithms  were  first  used  the  free 
surface  was  reconstructed  using  either  vertical  or  horizontal 
surfaces  [8,  3].  The  VOF  values  of  all  the  elements  are 
updated  each  time  step  by  calculating  the  mass  flux  between 
elements.  The  mass  flux  for  free-surface  elements  is 
calculated  by  taking  into  account  the  smaller  surface 
through  which  the  fluid  can  move  due  to  the  presence  of  the 
free  surface.  The  fixed  Cartesian  mesh  with  VOF  technique 
has  the  advantage  that  mesh  generation  is  straightforward. 
Also,  the  technique  can  easily  deal  with  floating  objects  and 
tank  baffles.  However,  the  disadvantage  of  the  technique  is 
that  the  fluid-solid  impenetrability  and  no-slip  boundary 
conditions  are  satisfied  only  in  a  time  average  sense.  Also, 
the  method  has  stability  and  accuracy  problems  when  the 
cut-cell  elements  at  the  solid-fluid  interface  become  small. 

(b)  Moving  ALE  mesh.  The  fluid  is  modeled  using  a  fluid  mesh 
that  moves  and  deforms  with  the  tank  as  well  as  with  the 
fluid’s  free-surface  using  the  ALE  formulation  [4,  6,  10]. 
The  main  disadvantage  of  this  method  is  that  it  does  not 
allow  large  surface  deformation  including  surface  break-up 
and  merging  unless  the  fluid  domain  is  re-meshed  [11,  12]. 
If  re-meshing  is  used  frequently,  then  there  is  the  problem 
of  degraded  solution  accuracy  due  to  re-interpolation  of  the 
solution  field  onto  the  new  mesh. 
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(c)  Fixed-fluid  mesh  with  the  Navier-Stokes  equations  written 
in  a  reference  frame  fixed  to  the  tank  along  with  a  VOF 
free-surface  model  [4,  13]. 

(d)  Moving  ALE  mesh  with  a  VOF  free-surface  model.  The 
fluid  domain  is  modeled  using  a  mesh  that  moves  and 
deforms  with  the  tank.  Then  the  VOF  technique  outlined 
above  is  used  to  model  the  liquid  free-surface.  This  method 
does  not  suffer  from  the  cut-cell  solid-fluid  interface 
boundary-condition  problem  like  method  (a).  It  also  allows 
modeling  surface  break-up  and  merging  without  the  need 
for  re-meshing.  This  is  the  method  that  is  used  in  the 
present  paper. 

(e)  Moving  ALE  mesh  with  a  level-set  free-surface  model  [5], 
This  method  is  the  same  as  method  (d)  except  that  a  level- 
set  method  is  used  to  model  the  liquid  free  surface.  The 
method  uses  a  smooth  scalar  function  defined  at  every  point 
in  the  fluid  domain  which  specifies  the  signed  smallest 
Eucledian  distance  between  the  point  and  the  interface.  The 
evolution  of  the  scalar  function  is  governed  by  a  convection 
transport  equation  where  the  interface  is  moved  with  the 
fluid  velocity. 

(f)  Lagrangian  Particle  methods.  Lagrangian  particles  are  used 
to  model  the  fluid.  A  contact  model  between  the  fluid 
particles  and  the  tank  is  used  to  model  the  fluid-structure 
interaction.  The  main  disadvantage  of  those  types  of 
methods  is  the  large  number  of  particles  and  computing 
time  needed  to  accurately  predict  the  fluid  motion.  A 
special  type  of  this  class  of  methods  which  has  been 
successfully  applied  to  fluid-structure  interaction  problems 
is  the  particle  finite  element  method  (PFEM)  [6]  in  which 
the  particles  are  used  to  generate  a  polyhedral  finite  element 
mesh  every  time  step  using  an  extended  Delaunay 
tesselation.  The  solution  of  the  incompressible  Navier- 
Stokes  equations  is  then  carried  using  that  mesh.  The  PFEM 
requires  less  fluid  particles,  however,  the  tessellation  step  is 
computationally  intensive.  Particle  methods  can  easily 
handle  surface  break-up  and  merging,  floating  bodies,  and 
fluid-solid  impenetrability  boundary  condition. 

A  review  of  multibody  dynamics  modeling  techniques 
including  deformation  reference  frames,  treatment  of  large 
rotations,  discretization  techniques,  finite  elements,  constraint 
and  contact  modeling,  and  solution  techniques  is  presented  in 
[14],  Flexible  multibody  systems  involving  liquid  filled  tanks 
are  generally  ground,  marine,  air,  or  space  vehicles  where  the 
following  need  to  be  accurately  modeled: 

•  Large  arbitrary  rigid  body  translation  and  rotation. 

•  Rigid  bodies. 

•  Flexible  bodies  which  can  be  composed  of  beam,  shell  or 
solid  elements.  Those  include  for  example  a  shell  model  for 
a  flexible  tank. 

•  Joints  including  spherical,  revolute,  prismatic,  and  universal 
joints.  Joint  friction  and  clearances  also  need  to  be  modeled. 

•  Frictional  contact.  For  example  for  ground  vehicle 
applications  the  rolling  frictional  contact  of  tires  need  to  be 
accurately  modeled. 

•  Transmission  components,  clutches,  and  brakes  for  ground 
vehicle  applications.  All  those  components  involve  friction. 


In  order  to  satisfy  the  above  requirements  a  flexible  multibody 
dynamics  modeling  technique  with  the  following  characteristics 
is  used: 

•  An  explicit  time-integration  solver  accurate  for  long 
simulation  times  [15]. 

•  Total  Lagrangian,  total  displacement  equations  of  motion 
formulation  with  the  degrees  of  freedom  referred  to  a 
global  inertial  reference  frame  [15,  16,  17,  18], 

•  A  library  of  truss,  beam,  and  solid  nonlinear  finite 
elements  with  Cartesian  coordinate  degrees  of  freedom 
allowing  arbitrarily  large  element  rotations.  Those  include: 
o  Torsional-spring  type  3-node  beam  elements  [15,  16]. 
o  Natural-modes  eight-node  brick  elements  [17,  18]. 

Those  elements  can  also  be  used  to  model  shells  and 
beams.  One  element  through  the  thickness  is  sufficient 
to  accurately  model  the  membrane,  shear,  and  bending 
characteristics.  They  do  not  exhibit  locking  or 
spurious  modes  (widely  used  techniques  to  alleviate 
locking  such  as  hourglass  control  lead  to  elements  that 
do  not  maintain  solution  accuracy  over  very  long 
solution  times).  They  are  computationally  efficient. 
Assumed  strain  elements  of  comparable  accuracy  are 
more  computationally  expensive.  Any  material  law 
can  be  used  with  those  elements  including:  linear 
elastic,  hyper-elastic,  and  non-linear  laws. 

•  Penalty  formulation  for  modeling  joints  (spherical, 
re  volute,  cylindrical  and  prismatic)  [19]. 

•  Rigid  body  rotational  equations  of  motion  are  written  in  a 
body  (material)  frame,  with  the  resulting  incremental 
rotations  added  to  the  total  body  rotation  matrix  [19]. 

•  Normal  contact  modeled  using  a  penalty  formulation  [20, 
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•  Frictional  contact  modeled  using  an  accurate  and  efficient 
asperity-based  friction  model  [22]. 

•  Special  elements  for  modeling  wheels/pulleys  [20], 
sprockets  [23],  and  clutches  [24]. 

•  General  tire  model  [21],  The  model  includes  the  details  of 
the  tire  construction.  The  tire  rubber  is  modeled  using  brick 
elements.  The  bead,  tread  and  ply  are  modeled  using  beam 
elements  along  the  tire  circumference  and  meridian 
directions  with  appropriate  stiffness  and  damping 
properties.  Normal  contact  between  the  tire  and  the  wheel 
and  between  the  tire  and  the  pavement  is  modeled  using 
the  penalty  technique.  Friction  is  modeled  using  an 
asperity-based  approximate  Coulomb  friction  model.  The 
tire  inflation  pressure  is  modeled  by  applying  a  force 
normal  to  the  inner  surface  elements  of  the  tire  with  the 
out-of-equilibrium  force  and  moment  applied  to  the  wheel 
to  guarantee  self-equilibrium  of  the  tire  and  wheel  under 
the  pressure  load  [21]. 

•  General  contact  search  algorithm  that  find  the  contact 
penetration  between  finite  elements  and  other  elements  as 
well  as  general  triangle  and  quadrilateral  surfaces. 

Two-way  coupling  between  the  multibody  system  (vehicle) 
motion  and  the  fluid  is  achieved  by  satisfying  the  following 
conditions  at  the  solid-fluid  interface: 

•  The  fluid  velocity  normal  to  the  solid’s  surface  must  be 
equal  to  the  normal  solid  velocity. 
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•  The  fluid  velocity  tangent  to  the  solid  surface  can  range 
from  being  equal  to  the  tangential  velocity  of  the  solid 
surface  (no  slip  condition)  to  being  free. 

•  No  additional  energy  or  momentum  to  the  system  should 
be  introduced  at  the  interface. 

The  proper  boundary  condition  that  does  not  violate  the  above 
conditions  arises  from  the  use  of  Newton’s  equations  of  motion 
to  find  a  common  normal  acceleration  for  the  fluid  and  the  solid 
at  the  interface.  The  tangential  fluid  and  solid  accelerations  can 
range  from  being  the  same  (no-slip  condition)  to  being 
completely  decoupled. 

Finally  in  the  present  paper,  a  single  computational  code  which 
uses  a  time-accurate  explicit  solution  procedure  is  used  to  solve 
both  the  solid  and  fluid  equations  of  motion.  Many  commercial 
software  and  studies  on  modeling  liquid  sloshing  coupled  with 
solid  body  motion  use  two  codes  which  pass  the  interface 
forces  and  motion  back  and  forth  and  iterate  on  the  two  codes 
until  equilibrium  is  achieved  [e.g.  3,  13].  This  approach  adds 
additional  computational  burden  and  does  not  achieve  the  same 
accuracy  as  the  single  integrated  code  solution. 

The  rest  of  this  paper  is  organized  as  follows.  In  Sections  2  and 
3  the  equations  of  motion  for  the  solid  and  fluid  are  presented. 
In  Section  4  the  fluid-structure  coupling  model  is  presented.  In 
Section  5  the  VOF  free-surface  model  is  presented.  In  Section  6 
the  overall  explicit  solution  procedure  is  outlined.  In  Section  7 
numerical  simulations  of  a  multibody  test  fixture  with  a  liquid 
filled  tank  are  presented.  Finally,  in  Section  8  some  concluding 
remarks  are  offered. 

2.  SOLID  EQUATIONS  OF  MOTION 

In  the  subsequent  equations  the  following  conventions  will  be 
used: 

•  The  indicial  notation  is  used. 

•  The  Einstein  summation  convention  is  used  for  repeated 
subscript  indices  unless  otherwise  noted. 

•  Upper  case  subscript  indices  denote  node  numbers. 

•  Lower  case  subscript  indices  denote  vector  component 
number. 

•  The  superscript  denotes  time. 

•  A  superposed  dot  denotes  a  time  derivative. 

The  translational  equations  of  motion  are  written  with  respect 
to  the  global  inertial  reference  frame  and  are  obtained  by 
assembling  the  element  equations.  The  finite  elements  used 
here  use  only  translational  DOFs  with  no  rotational  DOFs.  This 
is  advantageous  in  terms  of  computational  efficiency,  accuracy, 
and  robustness  [14],  Those  equation  also  include  the  rigid-body 
(such  as  the  wheel)  translational  DOFs.  The  equations  can  be 
written  as: 

M,x',  =  F4+F’4  (1) 

where  t  is  the  running  time,  K  is  the  global  node  number  (no 
summation  over  K\  K=\ ~>N  where  N  is  the  total  number  of 
nodes),  i  is  the  coordinate  number  (;=1,2,3),  a  superposed  dot 
indicates  a  time  derivative,  MK  is  the  lumped  mass  of  node  K ,  x 
is  the  vector  of  nodal  Cartesian  coordinates  with  respect  to  the 
global  inertial  reference  frame,  and  x  is  the  vector  of  nodal 
accelerations  with  respect  to  the  global  inertial  reference  frame. 


Fs  is  the  vector  of  internal  structural  forces,  and  Fa  is  the  vector 
of  externally  applied  forces,  which  include  surface  forces  and 
body  forces. 

For  each  rigid  body,  a  body-fixed  material  frame  is  defined. 
The  rigid  body  is  represented  by  one  node  located  at  the  body’s 
center  of  mass,  which  is  also  the  origin  of  this  frame.  The  mass 
of  the  body  is  concentrated  at  the  node  and  the  inertia  of  the 
body  given  by  the  inertia  tensor  Iy  defined  with  respect  to  the 
body  frame.  The  orientation  of  the  body-frame  is  given  by  R'° 
which  is  the  rotation  matrix  relative  to  the  global  inertial  frame 
at  time  t0.  The  rotational  equations  of  motions  are  written  for 
each  rigid  body  with  respect  to  its  body-fixed  material  frames 
as: 

IJK,  =  T,'„  +  714  -  (4  x  UkA  )L  (2) 

where  IK  is  the  inertia  tensor  of  rigid  body  K,  $K.  and  Ohj  are 

the  angular  acceleration  and  velocity  vectors’  components  for 
rigid  body  K  relative  to  it’s  material  frame  in  direction  j,  TsKi 
are  the  components  of  the  vector  of  internal  torque  at  node  K  in 
direction  ;,  and  TaKi  are  the  components  of  the  vector  of  applied 
torque.  The  summation  convention  is  used  only  for  the  lower 
case  indices  /  and  j. 

The  trapezoidal  rule  is  used  as  the  time  integration  formula  for 
solving  equations  (1)  for  the  global  nodal  positions  x: 

x'k,  =  x'k?  +  0.5  At  (x'KI  +  x'Kf )  (3a) 

x'kj  =  44  +  05  A  t  (x'Kj  +  x'4' )  (3b) 

where  At  is  the  time  step.  The  trapezoidal  rule  is  also  used  as 
the  time  integration  formula  for  the  nodal  rotation  increments: 

4  =^-*+0.5  A  t(e'KJ+d'Kn  (4a) 

A01,,  =  0.5  At  (%,+%,»)  (4b) 

where  A0Kj  are  the  incremental  rotation  angles  around  the  three 
body  axes  for  body  K.  The  rotation  matrix  of  body  K  (RK)  is 
then  evaluated  using: 

R'K  =  R'-\R(A0'Ki)  (5) 

where  R(A0'Ki)  is  the  rotation  matrix  corresponding  to  the 
incremental  rotation  angles  from  Equation  (4b). 

The  explicit  solution  procedure  used  for  solving  equations  (1-5) 
along  with  constraint  equations  is  presented  in  Section  6.  The 
constraint  equations  are  generally  algebraic  equations,  which 
describe  the  position  or  velocity  of  some  of  the  nodes.  They 
include: 

•  Prescribed  motion  constraints: 

/({*},<)  =  0  (6) 

•  Joint  constraints: 

/(M)  =  0  (7) 

•  Contact/impact  constraints: 

/((*})  >0  (8) 

The  penalty  technique  is  used  for  imposing  the  constraints  in 
which  a  normal  reaction  force  is  generated  when  a  node  is 
penetrates  in  a  contact  body  whose  magnitude  is  proportional  to 
the  penetration  distance  [20,  21].  An  asperity-spring  friction 
model  is  used  to  model  joint  and  contact  friction  [22]  in  which 
friction  is  modeled  using  a  piece-wise  linear  velocity-dependent 
approximate  Coulomb  friction  element  in  parallel  with  a 
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variable  anchor  point  spring.  The  model  approximates  asperity 
friction  where  friction  forces  between  two  rough  surfaces  in 
contact  arises  due  to  the  interaction  of  the  surface  asperities. 

3.  SEMI-DISCRETE  FLUID  EQUATIONS  OF  MOTION 

The  dynamic  response  of  the  fluid  is  described  by  the  ALE 
version  of  the  incompressible  Navier-Stokes  equations,  namely, 
the  equations  of  conservation  of  momentum  and  mass  for  a 


moving  deforming  control  volume:  which  are: 

\p*-  dV  =  f -u,  *dV  +  f'  ^  ^  V  + 

jpfdV 

V 

&dv+  f^-)dv  =  Q 

•  At  •  Ac, 

(10) 

P  =  Pa  +  rP 

(11) 

«i=«i  ~vi 

(12) 

TiJ=ADki5iJ+2(p  +  pl  +  p*)D.J 

(13a) 

Dtj  =  0.5{ASj/Acj  +  Sii  Jckt ) 

(13b) 

where  V  is  the  element  volume,  t  is  the  running  time,  p  is  the 
density  of  the  fluid,  u  is  the  fluid  velocity  vector  relative  to  the 
global  reference  frame,  u  is  the  fluid  velocity  vector  relative  to 
the  moving  fluid  mesh,  v  is  the  velocity  of  the  fluid  mesh,  P  is 
the  relative  pressure,  x  is  the  position  vector,  r  is  the  deviatoric 
stress  tensor,  D  is  the  rate  of  deformation  tensor,  f  is  the  body 

force  vector,  r  is  the  artificial  compressibility  parameter,  pn  is 
the  nominal  fluid  density,  p  is  the  fluid  viscosity,  p,  is  an 
additional  turbulence  viscosity  (set  proportional  to  vorticity 
magnitude),  and  p.  is  an  additional  free-surface  viscosity  used 
to  model  the  interaction  between  the  liquid  and  the  gas  that  is 
above  it.  Incompressible  flow  is  modeled  using  the  artificial 
compressibility  technique  [25].  A  finite  element  formulation  is 
used  to  derive  the  element’s  semi-discrete  equations  of  motion 
from  the  governing  equations  (9-13).  8-node  hexahedral 
elements  are  used  with  tri-linear  equal-order  velocity  and 
pressure  interpolation.  A  pressure  averaging  algorithm  [26]  is 
used  to  eliminate  pressure  checker-boarding  (due  to  the  use  of 
an  equal  order  interpolation  for  pressure  and  velocity).  The 
element  equations  are  assembled  into  the  global  semi-discrete 
equations  of  motion: 

MfKu'KI=F/Ki  (14) 

VfKK=QfK,  <15) 

where  Mfk  is  the  lumped  fluid  mass  of  node  K,  u'Kj  is 

component  i  of  the  fluid  acceleration  at  node  K,  F/'a  is 
component  i  of  the  fluid  forces  at  node  K,  Vkk  is  the  lumped 
fluid  volume  at  node  K,  P'  is  component  /  of  the  fluid  pressure 

rate  at  node  K  and  Qf'K.  is  component  i  of  the  fluid  pressure 

fluxes  at  node  K.  Those  equations  are  integrated  using  the 
trapezoidal  rule  along  with  an  explicit  solution  procedure  to 
yield  the  nodal  fluid  velocity  and  pressure: 

+0.5At(u'Ki+u'KiAl)  (16) 

^=V+0.5M4+^")  (17> 


4.  FLUID-STRUCTURE  COUPLING  MODEL 

Newton’s  equations  of  motion  are  used  to  find  a  common 
normal  acceleration  for  the  fluid  and  the  solid  at  the  interface. 
This  is  done  for  each  node  at  the  fluid-structure  boundary  as 

follows: 

(ms  +  m/)u„  —  Fluid  Forces + ^  StructureForces  (18a) 

(ms  +  mj)\>n  =  Fluid  Forces  +  StructureForces  (18b) 

where  ms  is  the  solid  mass  of  the  node,  mf  is  the  fluid  mass  of 
the  node,  u„  and  v„  are  respectively  the  fluid  and  solid 
accelerations  of  the  node  normal  to  the  fluid-structure  interface. 
The  tangential  fluid  and  solid  accelerations  (u,,v,)  are 
calculated  using  the  following  equations: 

((1  -  5)  ms  +  mr) iu  =  (1  -  StructureForces + ^  Fluid  Forces  (19a) 

(m,  +  (I  -  s)m,)i’,  =  V  StructureForces  +  (1  -  s)^  Fluid  Forces  (19b) 
where  s  is  the  slip  factor.  A  no-slip  condition  corresponds  to  a 
slip  factor  of  zero.  The  slip  factor  determines  how  much  of  the 
fluid  and  structure  forces  are  mutually  exchanged.  Equations  18 
and  19  are  written  for  all  fluid-solid  interface  nodes.  The  fluid 
mesh  must  move  and  deform  with  the  tank.  This  is  done  by 
modeling  the  fluid  mesh  using  very  light  and  compliant  (3 
orders  of  magnitude  less  than  the  tank)  solid  brick  elements 
(called  “mock”  mesh).  The  ALE  formulation  is  used  to  account 
for  the  fluid  mesh  deformation/motion. 


5.  VOF  FREE-SURFACE  MODEL 


For  each  fluid  element  a  VOF  value  between  0  and  1  is  defined, 
where  0  corresponds  to  empty  elements  and  1  corresponds  to 
elements  completely  filled  with  fluid.  The  elements’  VOF 
values  are  updated  each  time-step  by  moving  fluid  from  a 
completely  or  partially  filled  “donor”  element  to  an  empty  or 
partially  filled  neighboring  “acceptor”  element  using  the 
following  model: 

V,„  =  V-  VOF,  (20) 

Vm  =  V„(]~  VOF„)  (21) 


AtSAn.u  Vso  >  AV  and  Vm>Vm 


AV  = 


Vo 

Vm 


V-o  <AV 
Vm  <  AV 


(22) 


where  Vc  is  the  volume  of  the  element;  Vn  is  the  volume  of  the 
neighboring  element;  Vm  is  the  volume  of  the  element  occupied 
by  the  fluid;  Vm  is  the  volume  of  the  neighboring  element 
available  to  receive  fluid;  AV  is  the  volume  flow  through  the 
boundary  between  the  two  elements  in  a  time  step;  At  is  the 
solution  time  step;  S  is  the  surface  area  between  the  two 
elements;  A  is  a  value  between  0  and  1  indicating  the  free- 
surface  aperture  through  which  the  fluid  can  move  from  the 
element  to  the  neighboring  element;  n  is  the  unit  normal  to  .S'; 
and  u  is  the  fluid  velocity  vector  at  the  surface  .S’.  If  AV  is  less 
than  0  then  the  element  is  an  acceptor  element  and  the  VOF 
values  are  not  updated  because  they  will  be  updated  later  when 
the  neighbor  element  is  set  to  be  the  donor  element.  So  if  AV  is 
greater  than  0  then  the  VOF  values  are  updated  using  the 
following  equations: 

VOF,  =  VOF,  -  AT  IVc  (23) 

VOF„  =  VOF„  +  A  VI  V„  (24) 

The  free-surface  apertures  A  at  the  element  interfaces  are  used 
to  limit  the  fluid  flow  based  on  the  location  of  the  free  surface 
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inside  the  element.  A  is  calculated  as  follows.  If  the  VOF  value 
of  the  element  is  1  then  there  is  no  free-surface  at  the  element, 
therefore  A- 1.  For  elements  with  a  VOF  value  less  than  1,  the 
following  steps  are  used  to  calculate  A: 

•  Calculate  the  normal  to  the  surface  by  looking  at  a  stencil  of 
neighboring  elements  around  the  element.  This  is  done 
using  the  following  equation: 

nci=VO¥„kS„kn„ki  i=  1,2,3  (25) 

where  >ie,  is  ;th  component  of  the  normal  to  the  free-surface 
at  the  element,  VOF»*  is  the  VOF  value  for  neighboring 
element  number  k,  is  the  area  of  the  intersection  surface 
between  the  element  and  neighboring  element  k ,  and  n„k,  is 
the  component  i  of  the  normal  to  the  surface  between  the 
element  and  neighboring  element  k.  He  is  then  normalized 
into  a  unit  vector.  Figure  1  shows  a  2D  4-node  quadrilateral 
and  the  free-surface  along  with  the  normal  He . 

•  Calculate  the  apertures  A  for  each  neighboring  element  by 
constructing  a  planar  surface  with  normal  He  and  with  total 
volume  equal  to  VOF.  Ve  (see  Figure  1). 


Figure  1  Stencil  of  neighboring  elements  used  to  determine  the 
free-surface  normal  He  and  the  liquid  free-surface. 

6.  EXPLICIT  SOLUTION  PROCEDURE 

The  solution  fields  for  modeling  the  solid,  fluid  and  liquid  free- 
surface  are  defined  at  the  model  nodes.  These  are: 

•  Solid  translational  positions. 

•  Solid  translational  velocities. 

•  Solid  translational  accelerations. 

•  Solid  rotation  matrices. 

•  Solid  rotational  velocities. 

•  Solid  rotational  accelerations. 

•  Fluid  velocities. 

•  Fluid  accelerations. 

•  Fluid  pressure. 

•  Fluid  pressure  rate. 

•  Volume-of-fluid. 

•  In  addition,  eddy  kinetic  energy  and  eddy  kinetic  energy 
rate  can  be  added  to  model  turbulence  if  an  LES  turbulence 
model  is  used. 

The  explicit  time  integration  solution  procedure  for  modeling 
the  coupled  response  of  the  solid  (multibody  system),  fluid,  and 


liquid  free-surface  (using  the  VOF  formulation)  predicts  the 
time  evolution  of  the  above  response  quantities.  The  procedure 
is  implemented  in  the  DIS  [27]  (Dynamic  Interactions 
Simulator)  commercial  software  code  and  is  outlined  below: 

1-  Prepare  the  run: 

a.  Set  the  initial  conditions  for  the  solution  fields  identified 
above. 

b.  Create  a  list  of  all  the  finite  elements  (including  both  solid 
and  fluid  elements). 

c.  Create  a  list  of  elements  that  will  run  on  each  processor. 
This  is  done  using  an  algorithm  which  tries  to  make  the 
computational  cost  on  each  processor  equal. 

d.  Create  a  list  of  all  the  constraints  (including  both  solid  and 
fluid  constraints). 

e.  Calculate  the  solid  masses  for  each  finite  element  node  by 
looping  through  the  list  of  finite  elements.  Note  that  the 
solid  masses  are  fixed  in  time. 

f.  For  each  node  create  a  list  of  comer  and  edge  nodes  that 
are  connected  to  it  using  fluid  volume  elements. 

g.  VOF  preparations: 

i.  Find  a  list  of  the  volume  fluid  elements. 

ii.  Create  a  list  of  fluid  volume  elements  that  will  run  on 
each  processor.  This  is  done  using  an  algorithm  which 
tries  to  make  the  computational  cost  on  each  processor 
equal. 

iii.  For  each  element  find  all  neighboring  elements. 

iv.  For  each  element  find  the  element  VOF  using  the 
nodal  VOF. 

v.  Re-interpolate  the  elements’  VOF  to  nodal  VOF. 

h.  Loop  over  all  the  elements  and  find  the  minimum  time 
step  for  the  explicit  solution  procedure. 

i.  Loop  over  all  the  elements  and  create  a  list  of  wall  nodes. 
For  each  wall  node  find  the  list  of  fluid  boundary 
elements. 

2-  Loop  over  the  solution  time  and  increment  the  time  by  At 

each  step  while  doing  the  following: 

a.  Set  the  nodal  values  at  the  last  time  step  to  be  equal  to  the 
current  nodal  values  for  all  solution  fields. 

b.  Do  2  iterations  (a  predictor  iteration  and  a  corrector 
iteration)  of  the  following: 

i.  Initialize  the  nodal  fluxes  to  zero.  Those  include:  solid 
forces,  solid  moments,  fluid  forces,  boundary  fluid 
forces,  and  pressure  fluxes.  In  addition,  the  lumped 
nodal  fluid  volume  and  fluid  mass  vectors  are  also 
initialized  to  zero. 

ii.  Calculate  the  nodal  solid  and  fluid  fluxes  and  the 
lumped  fluid  volume/mass  vectors  by  looping  through 
all  the  elements  while  calculating  and  assembling  the 
element  nodal  fluxes  and  vectors.  This  is  the  most 
computational  intensive  step.  This  step  is  done  in 
parallel  by  running  each  list  of  elements  identified  in 
step  1  .c  on  one  processor. 

iii.  Find  the  nodal  values  at  the  current  time  step  using  the 
semi-discrete  equations  of  motion  and  the  trapezoidal 
time  integration  rule  (Equations  1-5  and  14-17). 

iv.  Execute  the  solid  and  fluid  constraints.  The  constraints 
set  the  nodal  value(s)  to  prescribed  values. 

v.  Apply  fluid-structure  interface  boundary  conditions 
for  all  wall  nodes  found  in  Step  Li  (see  Equations  18, 
19).  This  is  done  by  doing  the  following  for  each  wall 
node: 
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1 .  Find  the  normal  to  the  surface  at  the  wall  node. 

2.  Normalize  the  surface  normal. 

3.  Find  the  solid,  fluid  bulk  and  fluid  boundary 
forces  in  the  directions  normal  and  tangent  to  the 
surface. 

4.  Find  the  normal  and  tangential  solid  and  fluid 
accelerations  using  the  trapezoidal  integration 
rule  and  the  wall  slip  percentage. 

vi.  Set  the  pressure  boundary  conditions  at  the  free 
surface. 

vii.  Update  the  VOF  field: 

1.  For  each  fluid  element  calculate  the  element 
volume.  This  step  is  done  in  parallel  using  the  list 
of  fluid  elements  for  each  processor  found  in 
Step  l.g.ii. 

2.  For  each  fluid  element  find  the  apertures  through 
which  the  fluid  convects  to  each  neighboring 
element.  This  step  is  done  in  parallel  using  the 
list  of  fluid  elements  for  each  processor  found  in 
Step  l.g.ii. 

3.  For  each  fluid  element  use  the  apertures,  the 
element  volume,  the  element  current  VOF  value, 
and  the  element  nodal  velocities  to  update  the 
VOF  value  of  all  neighboring  elements  by 
finding  the  volume  of  fluid  that  left  the  element 
during  that  time  step  using  Equations  20-24.  This 
step  is  done  in  parallel  using  the  list  of  fluid 
elements  for  each  processor  found  in  Step  l.g.ii. 
Note  that  this  step  depends  on  the  order  of  the 
elements  in  the  list  of  elements.  However,  since 
the  updates  of  the  VOF  field  between  solution 
time  steps  are  small,  therefore  this  dependence  is 
generally  very  small.  In  order  to  assure  minimum 
dependence  on  the  elements’  order,  at  a  time  step 
the  elements  are  updated  from  first  to  last,  then  at 
the  next  time  step  they  are  updated  from  last  to 
first. 

viii.  Average  the  fluid  pressure  (This  step  eliminates  the 
pressure  checker-boarding  effect  and  allows  use  of 
equal  order  interpolation  for  both  pressure  and 
velocity). 

ix.  Go  to  the  beginning  of  step  2. 

An  advantage  of  explicit  solution  procedures  is  that  they  are 
“embarrassingly”  parallel.  The  above  procedure  achieves  near 
linear  speed-up  with  the  number  of  processors. 

7.  NUMERICAL  SIMULATIONS 

Figure  2  shows  snapshots  of  liquid  sloshing  in  an  oval 
tank.  The  liquid  modeled  is  water  (p  =  1 000  Kg, 
p  =  0.001  Kg/(m.sec)).  The  tank  dimensions  are:  0.324  m 
length,  width  0.428  m  and  height  0.276  m.  The  fluid  is  modeled 
as  incompressible  using  the  artificial  compressibility  technique 
with  an  artificial  sound  speed  factor  of  0.1  (i.e.  the  artificial 
sound  speed  in  the  water  is  taken  as  1483m/sec  xO.l  =  148.3 
m/sec).  Due  to  the  use  of  large  elements  near  the  solid  surface, 
full  slip  boundary  condition  at  the  wall  is  used.  Thus,  the 
viscous  wall  friction  effects  are  assumed  to  be  negligible.  The 
liquid  starts  from  the  initial  conditions  at  time  0  shown  in 
Figure  2  and  then  sloshes  in  the  tank  under  a  gravity  field  in  the 
vertical  direction  of  9.8  m/sec2. 
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Figure  2  Snapshots  of  the  DIS  simulation  of  liquid  sloshing  in  an 
oval  tank  from  time  0  to  0.6  sec. 


Figure  3  Multibody  test  fixture  with  a  tank  mounted  on  a 
suspension  system  composed  of  3  linear  springs. 

Next  the  above  tank  is  mounted  on  a  suspension  system  which 
is  mounted  on  a  test  fixture.  The  test  fixture  was  physically 
constructed  to  validate  the  computational  code.  Physical 
experiments  to  validate  the  computational  code  presented  in 
this  paper  will  be  conducted  in  the  near  future.  The  test  fixture 
simulates  some  of  the  maneuvers  that  a  typical  tank  in  a  tanker 
truck  is  subjected  to  such  as  roll  and/or  pitch.  The  test  fixture  is 
composed  of  the  following  components  (Figure  3): 

•  Three  linear  springs.  Simulate  the  suspension  system  of  the 
tanker  truck. 
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A  chassis.  The  springs  are  mounted  on  the  chassis  using 
spherical  joints. 

Two  linear  actuators  connected  to  the  chassis  provide  the 
ability  to  move  the  chassis  in  order  to  simulate  roll,  pitch 
and  combined  roll/pitch. 

A  rigid  grounded  base. 

Connecting  rods.  Provide  lateral  and  longitudinal  stability 
for  the  tank. 

1 7  spherical  joints. 

5  prismatic  joints  located  at  the  3  suspension  system  springs 
and  the  two  actuators. 


•  The  tank  is  assumed  to  be  flexible  and  is  modeled  using 
shell  elements. 

The  same  initial  conditions  as  in  Figure  2  where  tank  is  half¬ 
full  and  the  liquid  free-surface  is  vertical  are  used  as  initial 
conditions  for  fully  coupled  multibody-liquid  sloshing 
simulation  of  the  test  fixture.  Figure  4  shows  snapshots  of  the 
motion  of  the  test  fixture  and  the  liquid  from  time  0  to  0.75  sec. 
Figure  5  shows  snapshots  from  different  angles  taken  at  time 
0.5  sec. 
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8.  CONCLUDING  REMARKS 

A  finite  element  model  for  predicting  the  fully  coupled 
dynamic  response  of  flexible  multibody  systems  and  liquid 
sloshing  in  containers  was  presented.  The  model  has  the 
following  characteristics: 

•  Parallel  explicit  time-integration  solver. 

•  Library  of  accurate  large  rotation  finite  elements  including: 
truss,  beam,  shell  and  solid  elements.  The  elements  only  use 
Cartesian  coordinates  as  DOFs. 

•  The  fluid  mesh  is  modeled  using  a  very  light  and  compliant 
solid  mesh  which  allows  the  fluid  mesh  to  move/deform 
along  with  the  tank  using  the  Arbitrary  Lagrangian-Eulerian 
formulation. 

•  Acceptor-donor  VOF  algorithm  for  modeling  the  fluid's 
free-surface. 

•  The  motion  of  the  solid  and  fluid  is  referred  to  a  global 
inertial  Cartesian  reference  frame. 

•  A  total  Lagrangian  deformation  description  is  used  for  the 
solid  elements. 

•  The  penalty  technique  is  used  to  model  the  joints. 

Numerical  simulations  were  presented  for  a  half-filled  tank 
supported  by  linear  springs  mounted  on  a  test  fixture.  Physical 
experiments  will  be  conducted  in  the  near  future  to  validate  the 
computational  model  presented  in  this  paper. 
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